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Rainfall triggers more deep-seated landslides than 
Cascadia earthquakes in the Oregon Coast Range, USA 


S. R. LaHusen"*, A. R. Duvall’, A. M. Booth’, A. Grant’, B. A. Mishkin’, D. R. Montgomery’, 


W. Struble®, J.J. Roering®, J. Wartman7 


The coastal Pacific Northwest USA hosts thousands of deep-seated landslides. Historic landslides have primarily 
been triggered by rainfall, but the region is also prone to large earthquakes on the 1100-km-long Cascadia Sub- 
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duction Zone megathrust. Little is known about the number of landslides triggered by these earthquakes because 
the last magnitude 9 rupture occurred in 1700 CE. Here, we map 9938 deep-seated bedrock landslides in the Oregon 
Coast Range and use surface roughness dating to estimate that past earthquakes triggered fewer than half of the 
landslides in the past 1000 years. We find landslide frequency increases with mean annual precipitation but not 
with modeled peak ground acceleration or proximity to the megathrust. Our results agree with findings about other 
recent subduction zone earthquakes where relatively few deep-seated landslides were mapped and suggest 
that despite proximity to the megathrust, most deep-seated landslides in the Oregon Coast Range were triggered 


by rainfall. 


INTRODUCTION 

Landslides are one of the greatest secondary hazards of large earth- 
quakes. Hillslope failures triggered by shaking can affect homes and 
infrastructure and dam rivers and cause flooding, sometimes in cat- 
astrophic outbursts (1-4). In the Pacific Northwest of the United 
States, magnitude 8 (M8) to M9 earthquakes occur every 300 to 
500 years along the Cascadia Subduction Zone (CSZ) (5), a 1100-km- 
long megathrust fault that accommodates subduction of the Juan de 
Fuca Plate beneath the North American plate (Fig. 1). Empirical 
estimates for the number of triggered landslides for earthquakes of 
this magnitude range from 50,000 to 1,000,000+, with total predicted 
volumes of 10 to 110 km’ (6, 7), representing a risk to the millions of 
people living in the region. However, predictive models for earthquake- 
induced landslides from a potential CSZ event have high uncertainty 
and often fail to match the observations from similar subduction 
zone events. The few complete coseismic landslide inventories in 
subduction zones challenge assumptions about the number and vol- 
ume of landslides triggered by these earthquakes. 

Recent work has improved earthquake-landslide volume scaling 
relationships by accounting for rupture depth and topographic 
steepness; these models predict total coseismic landslide volumes of 
0.6 to 2.0 km? for moment magnitude (M,,) 9.0 CSZ earthquakes 
with rupture depths of between 20 and 10 km, respectively (8). Even 
these improved models overestimate the total coseismic landslide 
volume for recent subduction zone earthquakes by ~2 orders of 
magnitude, which may be an effect of the large distance between 
rupture hypocenters and subaerial landslide-prone topography. Only 
1226 landslides were mapped following the 2010 M, 8.8 Maule 
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earthquake in Chile (9). The following year, the 2011 M,, 9.0 Tohoku 
earthquake struck Japan, causing peak ground accelerations (PGAs) 
over 1.0g at many sites but resulting in just 3477 mapped landslides 
(10). Most of the landslides triggered in both earthquakes were shal- 
low, disrupted slides (defined here as landslides where the mobilized 
slide material mostly or entirely disintegrates during failure, resulting 
in flow or avalanche-style behavior) and lateral spreads. In Cascadia, 
the timing of the last great CSZ earthquake is well constrained to 
1700 CE (11), yet no prehistoric landslides have been definitively 
dated to that time (12). Recent landslide dating efforts have precisely 
determined the age of several landslide-dammed lakes in western 
Oregon, but none date to 1700 CE (13). Complicating the search for 
past coseismic landslides in Cascadia is the characteristically wet 
climate of the Pacific Northwest, where heavy rainfall often triggers 
landslides (14-16). Recent work in the Nepalese Himalaya, where 
landslides are triggered by earthquakes and monsoonal rain, suggests 
that rainfall is responsible for between 40 and 90% of the total land- 
slide erosion on 10*-year time scales (17). Because of the unique cli- 
mate and tectonic setting of Cascadia, the challenge of parsing the 
relative importance of these triggers remains essential to understand- 
ing what controls landslide hazards and, on longer time scales, land- 
scape denudation in the region. 

One barrier to understanding the behavior and triggering mech- 
anisms of past landslides is the difficulty of dating landslides on a 
regional scale. Landslide mapping inventories, however accurate, are 
of limited use in their ability to assess landslide triggering or event 
frequency because they often contain no age information for indi- 
vidual slope failure events. Recent work in Washington state has 
shown that for landslides in glacial sediment, surface roughness 
measured from lidar data can be used as a proxy for landslide age, 
enabling regional landslide chronologies to be constructed from a 
limited set of absolute dates (18, 19). This technique takes advantage 
of the topographic smoothing of landslide deposits over time as 
erosion and soil transport move material from convex areas to con- 
cave areas. Thus far, this technique has been applied to Quaternary 
glacial deposits (18, 19), and here, we demonstrate its effectiveness 
in dating bedrock landslides. 
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Fig. 1. Tectonic setting and geography of study area. Black outline depicts the 
central Oregon Coast Range (OCR) study area, which is underlain by the Eocene 
Tyee Formation and overlying Elkton Formation. The CSZ runs roughly north to south 
for over 1100 km just offshore. Shaded relief base map compiled from U.S. Geological 
Survey 30-m elevation data and National Oceanic and Atmospheric Administration 
bathymetry by the state of Oregon. Depth contours show the interface between the 
Juan de Fuca Plate and the overlying North American plate (27, 28). 


In this study, we apply surface roughness dating to an inventory of 
9938 manually mapped, deep-seated bedrock landslides across 
15,000 km” in the sandstones and siltstones of the Tyee and Elkton 
Formations (20) of the central Oregon Coast Range (OCR) (Fig. 2). 
Landslide deposits were mapped from 0.9-m bare-earth lidar imagery 
based on the presence of an arcuate headscarp and adjacent hum- 
mocky topography (fig. S1). The OCR, a mountainous swath of 
western Cascadia, has a well-documented abundance of prehistoric 
and historic landslides and is prone to both deep-seated bedrock 
landslides and shallow translational slides and debris flows (21-26). 
Thousands of landslides have previously been identified in the OCR 
(26), and automated landslide identification techniques suggest that 
up to 25% of the Tyee Formation is affected by deep-seated slides, 
underscoring the importance of this process to landscape evolution 
and morphology (24). Valley-spanning prehistoric landslides exist 
throughout the region, with remnant portions of deposits causing 
substantial upstream aggradation (22). Historic deep-seated landslides 
include a failure that occurred after days of heavy rain in 1975, dam- 
ming Drift Creek and creating a lake that persists today (24). 
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Numerical earthquake models suggest that the OCR, which over- 
lies the CSZ and in some places is less than 100 km from the associated 
offshore trench, may experience PGAs >0.6g during M9 earthquakes 
(27, 28 Fig. 1). When these strong ground motions interact with steep 
slopes often exceeding 30° in our study area, many models would 
predict widespread coseismic landslide triggering [e.g., (29)]. How- 
ever, the lack of historic records of coseismic landslides associated 
with the CSZ inhibits our understanding of where and how slopes 
fail during these large earthquakes. With the recent availability of 
<1-m-resolution bare-earth lidar coverage of the entire OCR, far 
more accurate and complete landslide mapping is now possible, as 
are surface roughness dating techniques that require high-resolution 
topographic data. Using 14 bedrock landslides with independent age 
constraints ranging from 7 + 1 years to 41,800 + 1000 years before 
present (ybp; here 2019) from this and other studies, we find an 
exponential relationship between landslide deposit surface roughness 
and age. We then use that calibrated roughness-age function to esti- 
mate the ages of all mapped landslides in the study area. With this 
landslide chronology, we first test whether predicted strong ground 
motions from CSZ earthquakes or mean annual rainfall exert con- 
trol on the spatial distribution of deep-seated slope failures over the 
past 1000 years. We then examine temporal trends in landslide fre- 
quency over the same time span to determine the importance of past 
large earthquakes as landslide triggers. 


RESULTS 

Surface roughness dating of bedrock landslides 

Previous work has tested the ability of numerous roughness metrics 
to both track known landslide ages and correctly identify relative 
ages among landslides with cross-cutting relationships (18, 19). For 
this study, we adopt the measure with the highest combined perform- 
ance from these previous tests: the two-dimensional (2D) continuous 
wavelet transform (CWT) (19). We take the mean magnitude of the 
wavelet coefficient, which is equivalent to topographic curvature, 
over the entire mapped deposit for each of 14 landslides of known 
age to represent their roughness (table S1 and fig. S1). Similar to 
previous work in Washington state, the relationship between land- 
slide deposit roughness and age is well fit by an exponential decay 
function (18, 19). We fit a robust (least absolute error) exponential 
regression using these roughness-age data, which yields the follow- 
ing function: t = 1.428 x 10° x e~*!”°®, where t is the estimated land- 
slide age, in years before 2019, and R is the average magnitude of the 
wavelet coefficient at a 20-m length scale within a mapped landslide 
polygon, in units of m7! (Fig. 3). We use this function to estimate 
the ages and uncertainties of all 9938 landslides across the study area 
based on their measured roughness value (Fig. 2B). As a conservative 
estimate of uncertainty for any single estimated landslide age, we 
calculate prediction bounds at a 95% confidence level for all age 
values. Uncertainty in predicted ages increases with time. For land- 
slides with historic estimated ages (less than 200 ybp), 95% predic- 
tion bounds yield an age range of 0 to 1300 years. This age range 
grows to 0 to 2400 years for landslides estimated to be 1000 years 
old. For landslides estimated to be thousands to tens of thousands 
of years old, uncertainty is plus or minus a few thousand years. While 
this technique is therefore not appropriate for pinpointing the trig- 
gering time for an individual slope failure, it offers a way to examine 
broadscale landsliding patterns in space and time across a landscape 
when applied to large landslide inventories. 
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Fig. 2. Landslide locations within the study area. (A) Manually mapped deep-seated landslide deposits (n = 9938) from the study area (outline in Fig. 1) shown as black 
dots over a base map of elevation and shaded relief. Dated landslides shown as white circles with ID numbers (table $1). Inset map for (B) outlined in black. (B) Detailed 
map of landslide deposits and estimated age ranges. Deep-seated landslides shown as polygons draped over 0.9-m bare-earth lidar-derived hillshade, with estimated 


ages from deposit surface roughness identified by polygon color. 
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Fig. 3. Plot of landslide age versus deposit surface roughness. (A) Open circles 
represent ages of landslides with independent timing constraints from radiocarbon, 
dendrochronology, or repeat aerial imagery plotted against the surface roughness 
of those deposits, quantified by the average magnitude of the wavelet coefficient 
for each polygon using a CWT witha 20-m Mexican Hat wavelet. Exponential decay 
function fit to these data shown as red line, where tis the landslide age (years be- 
fore 2019) and Ris the landslide roughness (mean wavelet coefficient of mapped 
polygon). Functional bounds show 95% confidence of fitted regression, plotted as 
blue dashed lines. Prediction bounds show 95% confidence for landslide age pre- 
diction, plotted as solid black lines. (B) Same data points and regression plotted on 
a semilog axis to better illustrate variability in landslide ages for slides less than 
1000 years old. 


LaHusen et al., Sci. Adv. 2020; 6:eaba6790 16 September 2020 


Spatial relationships between rainfall, CSZ earthquakes, 
and landsliding 
To understand whether the spatial distribution of recent landslides 
is related to CSZ earthquakes or regional rainfall patterns, we iden- 
tify all landslides with estimated ages less than 1000 ybp. For each 
slide, we extract values for modeled PGA from simulations of an 
M9 CSZ event (27, 28), 3D distance to the potential CSZ megathrust 
rupture plane (27, 28), and mean annual precipitation [MAP; (30)] 
at the centroid point of each slide (Fig. 4). We then plot landslide 
frequency per unit area against each of these metrics and test the 
statistical significance of each correlation by examining linear regres- 
sions and their associated confidence intervals. This analysis offers 
insight into the relative importance of earthquakes versus rainfall in 
controlling spatial patterns in the frequency of deep-seated landslides. 
Although predicted PGA for a M,, 9.0 earthquake on the CSZ 
megathrust is substantial (ranging from 0.6g in the west to 0.2g in 
the east within the study area), our analysis shows no statistically 
significant linear correlation (Fig. 4A). A linear regression fit to the 
PGA-landslide frequency points yields the following equation with 
95% confidence intervals: Fis = 1.07 x 10-°(+1.169 x 10°) x PGA + 
5.79 x 10°>(+12.2 x 107°), where Fs is the yearly landslide frequency 
per square kilometer and PGA is unitless peak ground acceleration 
expressed as the proportion of g (R’ = 0.4). To account for the lack of 
any historic measurements of PGA from large-magnitude earthquakes 
along the CSZ, we also test these enigmatic results by examining a 
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Fig. 4. Plots showing the relationship between landslide frequency during the past 1000 years and PGA, distance to CSZ rupture surface, and mean annual 
rainfall. Area-normalized landslide frequency is plotted against (A) modeled peak ground acceleration (27, 28), (B) 3D distance to the modeled CSZ megathrust rupture 
surface (27, 28), and (C) MAP. Best-fit linear regression is shown as red line with shaded 95% functional confidence interval. Values for each of the three metrics considered 


were extracted at the centroid point of each landslide less than 1000 years old. 


simpler and more objective measure: 3D distance from each land- 
slide to the CSZ plate interface. We observe a similarly counterintui- 
tive pattern of landslide frequency weakly increasing with distance 
from the fault before dropping off at distances greater than 30 km, 
once again resulting in a lack of statistically significant linear trend: 
Frs = 9.64 x 10°°(+17.3 x 10°) x D + 1.48 x 10°°(+47.83 x 10°), 
where D = 3D distance from landslide center to the rupture plane, 
in kilometers (R’ = 0.2; Fig. 4B). Using 3D distance from a potential 
rupture surface accounts for the relatively shallow dip angle of the 
CSZ megathrust (27, 28), although a similar overall pattern is observed 
when 1D map distance to the offshore CSZ fault trace is considered. 
When mean annual rainfall is plotted against landslide frequency, 
we do find a statistically significant positive correlation in the follow- 
ing form: Fis = 12.3 x 10-°(+6.15 x 10°) x MAP + 5.79 x 10°-°(+12.2 x 
10°), where MAP is mean annual precipitation, in meters (R? = 0.7; 
Fig. 4C). The topography of OCR acts as an orographic barrier to 
offshore storms, creating a rain shadow that results in MAP that ranges 
from less than 1.5 m/year to over 2.5 m/year within the study area 
(30, 31). In wetter portions of the study area that receive ~2.5 m/year, 
landslide frequency is double that in the drier portions that receive less 
than 1.5 m/year (Fig. 4C). Collectively, these analyses give no indication 
that CSZ megathrust earthquakes control the location of landslides 
triggered during the past 1000 years but do show a strong correlation 
between rainfall and landslide frequency during the same period. 


Estimating landslide frequency through time in the 

central OCR 

The predominant pattern in landslide frequency within the study 
area is an apparent nonlinear decrease in frequency with age (Fig. 5A). 
This same pattern was observed in previous landslide chronologies 
(18, 19) and likely reflects both a real preservation bias of young slides 
and a bias effect associated with using an exponential function to 
estimate age (Supplementary Materials). The preservation bias should 
be expected as ongoing hillslope and fluvial processes tend to re- 
shape relict landslide deposits, and as old slides are reactivated. 
Equilibrium hillslope adjustment time scales calibrated to the OCR 
suggest that the topographic signature of the deepest landslides in 
the study area may persist for tens of thousands of years after failure 
(32), while the legacy of smaller, shallower landslides is more short 
lived. This is supported by the distribution of landslide deposit area 
for young versus old landslides: Deposits estimated to be tens of 
thousands of years old have nearly double the median and mean 
area compared with landslides in the past 1000 years (table S2). 
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Together, these processes lead to a preservation bias toward 
younger and larger landslides, similar to the well-documented ap- 
parent decline of sedimentation and erosion rates with geologic age 
that has been attributed to biased representation of the geologic re- 
cord (33). As a result, we cannot directly compare the true frequency 
of landsliding from tens of thousands of years ago to present day, 
but we can statistically identify peaks in the landslide age-frequency 
data that differ significantly from the background trend and may 
represent periods of widespread landslide triggering due to extreme 
storms or large earthquakes. To test whether there is a signal from 
past CSZ earthquakes in the distribution of landslide ages, we focus 
on the age-frequency data for the past 1000 years, where the onshore 
and offshore geologic records suggest that three full-length (~M8.7 
to 9.1) CSZ ruptures have occurred, including the most recent 1700 CE 
earthquake as well as earthquakes at 552 + 83 ybp and 868 + 58 ybp 
(5). After selecting a subset of 2676 landslides with estimated ages 
<1000 years, we bin landslide frequency into twenty 50-year intervals 
and then estimate slide frequency per year for each bin (Fig. 5B). 

Landslide frequency over the past 1000 years reveals a nonlinear in- 
crease in landslide frequency toward the present (to the left on Fig. 5), 
consistent with the overall trend observed for the past ~50,000 years. 
Superimposed on this trend are two sharp peaks at 125 and 225 ybp, 
and a broad peak centered at 375 ybp (Fig. 5B). To test whether these 
peaks represent noise in the dataset or a real increase in landslide 
frequency at these times, we compare this best-fit frequency plot with 
a suite of frequency plots generated from a bootstrap analysis. This 
bootstrap analysis allows the roughness and age values for landslides 
of known age to vary within their individual uncertainties and then 
fits a unique regression for each of 10* bootstrap runs. Given the 
limited number of landslides with age constraints used to fit our 
roughness-age function, this analysis helps to address error in our 
landslide frequency plots caused by uncertainty in the roughness-age 
regression. When we plot the bootstrapped mean frequency as well 
as the lower 5th and upper 95th percentiles, a trend of decreasing 
frequency with age is observed, similar to the frequency plot gener- 
ated by the best-fit regression; however, no peaks appear in the 
bootstrap mean distribution. When individual frequency plots from 
the bootstrap analysis are examined, peaks of similar amplitude as 
the peaks observed in the best-fit frequency plot appear at random 
times in the 1000-year history. Collectively, this analysis suggests 
that the sharp peaks observed in the best-fit frequency plot may not 
be statistically significant, and thus, deviations from a constant rate 
of landsliding during the past 1000 years should be viewed with 
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Fig. 5. Distribution of landslide ages and estimated landslide frequency. (A) Histogram of landslide ages for the past 50,000 years across the study area shows a 
nonlinear decrease with time, which is likely the result of the preservation bias of younger-landslide deposits. (B) Mean yearly landslide frequency estimated for the past 


1000 years. Frequency estimated from the best-fit roughness-age regression is shown as solid black line. A bootstrap analysis of 104 possible roughness-age regressions 
fit to our data is shown in blue, where the solid line is the mean and the dashed lines are the lower 5th and upper 95th percentiles of the bootstrap analysis. Timing of 


presumed full-rupture earthquakes shown as red solid lines with +o uncertainty (5). 


caution. If the peaks in the observed landslide frequency plot are 
noise, as the bootstrap analysis suggests, then there is no justification 
for inferring widespread landslide-triggering events at those times. 


Quantifying coseismic landslide triggering rates from 
CSZ earthquakes 
Although we do not detect clustered landslide activity in our analysis, 
the overall shape of the curve is useful in determining the relative 
importance of earthquakes as a triggering mechanism in the OCR. 
Because of the large uncertainty in the roughness-age technique, 
pulses of coseismic landslides may not show up as discrete peaks in 
the resulting landsliding frequency data but could form broad 
undulations that control the overall shape of the dataset over the 
past 1000 years. Moreover, there has not been a large earthquake in 
Cascadia for 319 years, so we may expect that the frequency plot 
from a landslide inventory triggered primarily by earthquakes would 
peak around the time of the 1700 CE earthquake and decrease to- 
ward present day. To systematically explore these scenarios and to 
place bounds on the number of possible earthquake-triggered land- 
slides in the past 1000 years, we generated simulated landslide 
chronologies and compared those with the observed landslide dataset. 
By superimposing pulses of coseismic landslides for each of the 
three large CSZ earthquakes in the past 1000 years (5) on a specified 
background rate of landsliding, we compared modeled frequency 
plots for three hypothetical landslide histories: (i) All landslides are 
seismically triggered; (ii) all landslides are triggered by background 
processes (we assume rainfall as the predominant background trigger); 
or (iii) an equal mix of seismic- and rainfall-triggered landslides. 
These simulated landslide chronologies allow us to discern the ex- 
pected form of landslide frequency plots for the past 1000 years for 
different triggering scenarios. By accounting for uncertainty in the 
roughness-age dating method that arises from natural variability in 
roughness values for landslides of the same age, we present a statis- 
tically justified way of interpreting our observed frequency data. 
The simulated landslide history composed entirely of “background” 
landslides, where no coseismic landslide pulses are added, most 
closely follows the observed frequency plot (Fig. 6A). This is expected 
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as the background rate and preservation bias in the simulated histories 
are tuned to match the total number of slides in the observed fre- 
quency data. For this “purely background” simulated history, the 
observed data fall almost entirely within the 5th and 95th percentile 
bounds of the model runs. As we increase the number of slides that 
are triggered in each earthquake pulse in the simulated landslide 
histories, the resulting frequency plots deviate further from the ob- 
served data. In model runs where 50% of the total landslides are 
coseismic, corresponding to 600 landslides triggered per earthquake, 
the observed data fall outside the modeled bounds for two distinct 
intervals in the 1000-year history: from ~700 to 250 ybp and from 
~150 ybp to present (Fig. 6B). The purely coseismic simulated history, 
for which 1300 landslides are triggered per earthquake, reveals the 
strongest deviation from the observed frequency data, with minimal 
overlap between simulated and observed landslide frequency (Fig. 6C). 
In the 100% coseismic simulation, the discrepancy between observed 
and simulated frequencies is caused by a shift in the overall shape of 
the curve, not by the formation of distinct frequency peaks at the 
times of earthquakes. The diffuse nature of these peaks causes them 
to overlap such that they are indiscernible. This effect is due to both 
natural variability in initial roughness conditions immediately after 
slope failure and variability in the quality of lidar data used to mea- 
sure roughness, which together result in a distribution of roughness 
values for landslides of the same age. When we shorten the band- 
width of this roughness distribution by decreasing the SD of rough- 
ness values used in the model runs, frequency peaks do appear (fig. S2). 
However, this requires unrealistic SD values, which are substantially 
lower than the measured SD and much lower than the more conserva- 
tive SD that we adopt for all model runs (see Materials and Methods 
and the Supplementary Materials). 


DISCUSSION 

In this study, we demonstrate that roughness dating of landslide 
deposits can be applied to large inventories of deep-seated bedrock 
landslides. In conjunction with previous studies of landslides in glacial 
sediments, this work suggests that roughness dating is a broadly 
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landslides mapped in the study area is plotted as a solid black line. Simulated landslide frequency for 104 modeled scenarios where no landslides are triggered during past 
full-rupture CSZ earthquakes is shown in red, where the solid line is the mean and the dashed lines are the lower 5th and upper 95th percentiles of all model runs. In the 
modeled scenarios, landslides are triggered at a steady rate. (B) Observed frequency shown as solid black line and simulated landslide frequency shown in red for modeled 
scenarios where 50% of the total modeled landslides are triggered by CSZ earthquakes, equivalent to 600 landslides per earthquake for three earthquakes during the past 
1000 years. The other 50% of landslides are triggered at a steady rate. Individual coseismic pulses associated with each earthquake are indiscernible due to roughness 
variability for landslides of the same age. (C) Observed frequency shown as solid black line and simulated landslide frequency shown in red for modeled scenarios where 


100% of the total modeled landslides are triggered by CSZ earthquakes, equivalent to 1300 landslides per earthquake. 


applicable tool suitable for use in diverse landscapes underlain by 
either bedrock or Quaternary sedimentary deposits (18, 19). In the 
Tyee Formation of the OCR, we show roughness dating to be a 
powerful tool for rapidly estimating the ages of nearly 10,000 land- 
slides across 15,000 km’, albeit with less accuracy than traditional 
landslide geochronometers such as radiocarbon dating. While this 
uncertainty may be too high to identify peaks in landslide frequency 
associated with individual earthquakes that occur every few hundred 
years in this region, the technique does offer an unprecedented means 
of parsing the relative importance of coseismic landslide triggering 
versus other triggering mechanisms. 

Because there has not been a large earthquake in the region in 
over 300 years, the shape of the observed versus modeled landslide 
frequency plots differs substantially, with frequency dropping off 
rapidly toward the present in simulated landslide histories that in- 
corporate a larger relative proportion of coseismic landslides. The 
observed landslide frequency over the past 1000 years shows the 
opposite: a slight increase in landslide frequency toward present day. 
Some of this increase is likely due to preservation bias, although we 
cannot discount the possibility of a real landslide frequency increase, 
potentially driven by widespread logging beginning in the late 19th 
century. The observed frequency is best reproduced with model 
simulations that include no pulses of coseismic landslides during 
the times of great CSZ earthquakes. In simulations where coseismic 
landslides account for half of the total landslides, the modeled fre- 
quency through time deviates substantially enough from the observed 
frequency that these simulations likely represent a conservative 
bound on the maximum proportion of coseismic triggering: equiv- 
alent to no more than 600 landslides per earthquake event. We cannot 
dismiss the possibility that CSZ earthquakes prime hillslopes for 
precipitation-induced landslides by damaging the uppermost few 
meters of regolith and bedrock without inducing catastrophic failure 
(17, 34). However, if this effect was substantial but short-lived, as is 
the case for earthquakes on shallow thrust faults (34), then we would 
still expect to see a peak in landslide frequency sometime after the 
last CSZ earthquake, but before present day. 

The apparent absence of coseismic landsliding in the study area 
may be unexpected, given the proximity of the region to the CSZ 
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megathrust and the intense ground shaking predicted by numerical 
earthquake simulation models (27, 28). However, these results are 
consistent with the conclusions drawn from recent coseismic land- 
slide inventories mapped after subduction zone earthquakes in Maule, 
Chile (9) and Tohoku, Japan (10). These studies found that the total 
number of coseismic landslides fell one to two orders of magnitude 
short of expected values based on earthquake magnitude-landslide 
frequency scaling relationships developed for crustal faults (7). More- 
over, coherent, deep-seated landslides made up the smallest number 
of landslides per failure style in the mapped inventories. Of the 1226 
landslides that Serey et al. (9) mapped after the Maule earthquake, 
only 8 were identified as coherent, deep-seated landslides. Improved 
predictive models that account for depth to the rupture surface still 
overpredict landslide volume for recent subduction zone earthquakes 
(8), although the relatively shallow dip angle of the CSZ megathrust 
positions a greater area of landslide-prone terrain directly above the 
potential rupture interface compared with other subduction zones. 
Still, we see no evidence that predicted PGA or 3D distance to the 
rupture interface exerts any control on spatial patterns in landslide 
frequency across the study area over the past 1000 years (Fig. 4). 
While catastrophic deep-seated landslides are large and may be 
very destructive, this study does not consider the potential for smaller 
shallow translational landslides or deep-seated slides that are only 
minimally displaced during a reactivation event. Topographic evi- 
dence for each of these landslide styles may not be perceptible in 
1-m lidar data, especially after hundreds of years, so we should con- 
sider the potential for these types of landslides to be triggered on a 
widespread scale during large earthquakes. Specifically, shallow 
landslides smaller than our threshold of 5 x 10° m? may be expected 
to occur in far greater number, on the basis of both landslide area- 
frequency scaling relationships (8, 35) and recent subduction zone 
earthquake-triggered landslide inventories (9, 10). Shallow slides and 
debris flows are particularly prevalent in steep catchments in the 
central OCR (36) and, importantly, could represent the dominant 
style of coseismic landslide in the study area. Colluvial hollows, a 
common debris flow source, refill with hillslope material, and shallow 
landslide deposits are less voluminous and more easily evacuated 
(37). This makes identifying and mapping these landforms difficult 
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or impossible hundreds of years after the last large earthquake. In 
contrast, scars and deposits associated with deep-seated landslides 
may last tens of thousands of years, a time span that may encompass 
dozens or even hundreds of earthquakes in seismically active regions. 

Given the well-established importance of rainfall as a landslide 
trigger (38-40), the lack of evidence for seismic triggering of deep- 
seated landslides, and the strong positive correlation between MAP 
and landslide frequency, we propose that precipitation likely triggers 
most landslides considered in this study. While we do not attempt to 
quantify erosion or rates of sediment delivery by rainfall-triggered 
landslides in the OCR, a recent study of the Nepalese Himalaya found 
that, if the 2015 Gorkha earthquake is representative of landslide- 
triggering earthquakes over 10°-year time scales, then earthquakes 
only trigger 10% of the total landslide volume in the region. If the 
entire suite of possible earthquake magnitudes and recurrence in- 
tervals are accounted for, then coseismic landslide volumes may 
account for up to 60% of the total. This attributes 40 to 90% of the 
total landslide erosional budget to rainfall-induced landslides, a 
similar conclusion to what we draw from the OCR, despite the sub- 
stantial differences in tectonic setting, lithology, and climate. Our 
conclusions are bolstered by recent efforts to pinpoint the timing of 
individual landslide deposits in the OCR using dendrochronology, 
often yielding ages with uncertainties of less than 1 year (13). This 
ongoing work has revealed multiple deep-seated landslides that date 
to 1890 CE, when historical records indicate severe flooding in the 
region (41), but has not identified any 1700 CE coseismic landslides. 
Together, both methods converge on a similar dearth of 1700 CE 
deep-seated landslides and suggest that the primary hazard posed 
by these types of landslides in Cascadia is driven by rainfall rather 
than large-magnitude subduction zone earthquakes. 


MATERIALS AND METHODS 

Landslide mapping and dating 

To construct a surface roughness—age relationship, we mapped 
9938 deep-seated landslides in the Tyee Formation and overlying 
Elkton Formation. Because surface roughness dating requires the 
assumption that landslides smooth over time because of erosion and 
soil diffusion, landslide deposit polygons had to be mapped carefully 
to avoid headscarps, incised stream gullies, and river cutbanks. Map- 
ping was completed manually using bare-earth lidar-derived hillshade 
and slope images. Only the presently intact surfaces of landslide de- 
posits were mapped to capture areas dominated by diffusive smooth- 
ing processes. Major roads and anthropogenic structures were removed 
from polygons before roughness was calculated (fig. $1). 

We mapped only coherent rotational and translational deep-seated 
landslide deposits with areas of at least 5 x 10° m” and with at least 
one dimension over 100 m. Although we can discern smaller land- 
slides in the lidar imagery, we use a 20-m wavelet to measure curva- 
ture for a roughness metric, so slides much smaller than this seemed 
inappropriate for consideration in this study. We define coherent 
landslides as slope failures where the mobilized material has not com- 
pletely disintegrated to the point of allowing flow or avalanche-style 
runout. Coherent slides more consistently preserve roughness ele- 
ments like intact blocks, minor scarps, hammocks, and compres- 
sional ridges needed for roughness dating. Wherever remobilized 
portions of older landslide deposits satisfying these criteria were ob- 
served, they were mapped as individual slides. To satisfy the assump- 
tion inherent to surface roughness dating that all landslides fall within 
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a range of expected initial roughness values directly after they occur, 
certain failure modes were excluded from the analysis. Landslides 
with channelized travel paths, debris flows, and earthflows were not 
considered, as the mechanics of these landslide types do not univer- 
sally produce classic deep-seated landslide morphologies such as 
hummocky topography, scarps, and displaced blocks characterized 
by high initial roughness, as is typical of coherent landslides. 

A suite of landslides with independently constrained ages is nec- 
essary for roughness-age dating, as these landslides provide the data 
to regress against roughness to define the roughness-age function 
used to estimate the ages of other landslides in the study area. We 
used radiocarbon dating and historical repeat aerial imagery to ascer- 
tain the timing of failure for six landslides in the mapped chronology 
(table S1). In addition, we used dates from previous work, including 
recent dendrochronology dating, for 8 other landslides (21, 22, 25, 26, 
41-43) yielding a total dataset of 14 landslides of known age (table 
S1). Because we date multiple landslides younger than 1950 CE, we 
report dated landslide ages in this study as years before 2019, as 
opposed to the radiocarbon convention of years before 1950. In this 
manuscript, we use the phrase “years before present” and its abbre- 
viation “ybp” to mean years before 2019. 


Surface roughness analysis and landslide age estimation 
For roughness dating to be applicable to a landslide-prone region, 
the initial roughness of landslide deposits is assumed to be similar, 
falling within a range of relatively high expected values. An additional 
assumption is that the material properties of the landslide and the 
climate of the study area are both similar enough that the rate of 
erosion and soil diffusion does not vary markedly among landslides 
within the inventory. Moreover, while this technique has been shown 
to be applicable to many landslides across large study areas (19), it 
has not yet been implemented on bedrock landslides. To test the 
utility of the approach on bedrock landslides, we focus our efforts on 
a swath of the central OCR, underlain entirely by Eocene marine sedi- 
mentary rock of the Tyee Formation and overlying Elkton Formation. 
We measured the roughness of each mapped landslide using the 
2D CWT with a Mexican Hat wavelet (fig. S1), which is equivalent 
to calculating topographic curvature (the Laplacian of elevation) at 
the specified length scale [(44) and the Supplementary Materials]. 
We first calculated the CWT at spatial scales of 4, 15, 20, and 45 m 
to find the scale that offered the best fit to the landslides of known 
age (table S3). For each CWT analysis, roughness values were com- 
puted for each cell in the 0.9-m bare-earth digital elevation model, 
and mean values of the magnitude for each landslide deposit of known 
age were extracted from mapped polygons and used to plot age against 
roughness. Two regression methods, robust (least absolute error or 
L1 norm) and least squares (L2 norm), were used to fit single exponen- 
tial functions to these roughness-age datasets. We tested the perform- 
ance of both methods in fitting the entire age range of dated landslides 
as well as a subset of young landslides (<1000 ybp). It was important to 
select a fit that did not result in a systematic overprediction or under- 
prediction in age estimations for the youngest landslides, as we used 
this fit to examine frequency trends in the past 1000 years. When the 
entire range of ages in the dataset was considered, the SE of the least 
squares regression was slightly lower (~1.2 x 10° years) compared 
with the robust regression (~1.6 x 10° years). The least squares re- 
gression optimized the fit for the oldest (>40,000 ybp) landslides but 
systematically overestimated the ages of the younger-landslide data 
points. When just the 10 landslides less than 1000 ybp were considered, 


7 of 10 


0202 ‘9¢ 42q0190 UO /bio BeWadUaIDS seouUeApe//:d}]4 WO’ papeojumMoq 


SCIENCE ADVANCES | RESEARCH ARTICLE 


the least squares regression resulted in nearly double the SE as the 
robust regression: ~400 years versus ~200 years, and examination 
of residual plots showed less bias in the robust regression. Because 
the robust regression method is less sensitive to outliers than the 
least squares regression, it most effectively fits a function across the 
age range of the dataset, especially for the younger (<1000 ybp) 
landslides. This robust regression, which we deem the best fit to our 
observed data, yields the following function that we use to estimate 
landslide age across the study area: t = 1.428 x 10° x e*!?°®, where 
t is the estimated landslide age, in years before 2019, and R is the 
average magnitude of the wavelet coefficient at a 20-m length scale 
within a mapped landslide polygon, m‘’ (Fig. 3). 


Analysis of PGA, distance to CSZ megathrust, and MAP 

To identify any correlation between the locations of recent landslides 
and patterns in modeled PGA, 3D distance to the CSZ, and MAP, 
we plot histograms of these data extracted from the center of each 
mapped landslide with estimated ages less than 1000 ybp. For mod- 
eled PGA, we adopt a mean of 30 simulated M9 CSZ earthquakes 
designed to capture a wide range of potential ground-shaking out- 
comes from variations in hypocenter location, down-dip rupture 
limit, slip distribution, and subevent locations (27, 28). This mean is 
downsampled to a 1-km grid using bilinear interpolation to smooth 
unrealistically large jumps between cells in the raw data. We also 
consider a simpler measure of CSZ seismic hazard: the minimum 
3D distance from each landslide point to the CSZ megathrust rupture 
surface. We take the rupture surface to be the interface between the 
Juan de Fuca Plate and the overlying North American plate from 
depths of 5 to 30 km (27, 28). For rainfall, we use a 30-year record of 
MAP (30). Although MAP does not capture the hourly intensity of 
individual storms, the deep-seated slides that we consider in this 
study are more likely to be sensitive to total precipitation volume on 
time scales of days to months [e.g., (39)]. Because of the highly sea- 
sonal nature of precipitation in the OCR along with the spatial vari- 
ability in winter storms, a 30-year record of MAP is a reasonable 
metric to capture the overall patterns in rainfall that are likely to 
control precipitation-induced landsliding. 

For each of the three metrics, we plot histograms by first dividing 
the middle 95% of the total range of values across the study area into 
10 equal intervals. Excluding the upper and lower 2.5% helps to 
avoid extreme values that affect little of the study area. The number 
of landslides for each value interval is then normalized by total area 
that interval represents within the study area and, lastly, converted 
into a yearly frequency. The resultant y axis shows the average fre- 
quency of landslides per kilometer for each equal interval of PGA, 
distance to CSZ, and MAP. To identify the correlation between 
each of the three metrics and the locations of recent landslides, we 
fit a line to the data in each plot and calculate the R’ value to test the 
goodness of fit (Fig. 4). 


Simulated coseismic landslide chronologies 

If large CSZ earthquakes do periodically trigger widespread land- 
sliding in the study area, then we would expect the resultant peaks 
in the histogram to be diffuse rather than sharp, because of the vari- 
ability of measured roughness values within each pulse of coseismic 
landslides. To account for this variability and to aid in a more robust 
interpretation of the observed landslide frequency data, we generated 
simulated landslide chronologies. These simulated chronologies offer 
a means to estimate the proportion of landslides in the study area 
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triggered by large CSZ earthquakes. Specifically, they allow us to 
better determine what a coseismic signal should look like. By compar- 
ing the observed and synthetic datasets, we can also deduce whether 
the apparent steady increase in landslide frequency that we observe 
could be due to uncertainty in the roughness dating technique. 

We assume that when a pulse of landslides occurs, such as from 
a large earthquake, landslides fail instantaneously, and their deposits 
have normally distributed roughness values. Over time, both the 
mean and the SD of the roughness values that compose this pulse 
decrease. The SD decreases because the roughest landslides have 
more regions of high curvature, and therefore smooth considerably 
faster than the landslides with lower initial roughness, effectively 
shrinking the distribution. This has been demonstrated numerically 
(19) and is supported by both geomorphic transport laws (45) and 
the exponential form of the roughness-age regression itself. To illus- 
trate this effect in the bedrock setting of the central OCR, we model 
the topographic evolution of the roughest and smoothest modern 
landslides in our dataset that have failed in the past 20 years. This 
exercise illustrates how the initial difference in roughness values be- 
tween the two slides rapidly shrinks after only a few hundred years 
(figs. S3 to $5). 

To calculate the SD of roughness values for a specific age, a dataset 
of landslides of known age would ideally include many landslides of 
the same age. From our landslides with known ages, we calculate 
the SD for four landslides with ages that fall within 100 years of the 
1700 CE earthquake, yielding a value of 8.4 x 10~* m™'. However, to 
address the uncertainty in calculating the SD from a small sample, 
we chose to double the measured SD and instead use an SD of 1.7 x 
10° for all modeled landslide distributions. Rather than attempt to 
scale down SD with age, we take a conservative approach to this 
uncertainty and use the same value for the entire 1000-year simula- 
tion. We acknowledge that the modeled roughness SD has a substan- 
tial effect on the resulting landslide frequency plots. Very high SD 
values yield estimated ages that are skewed toward present day, pro- 
ducing a steep nonlinear increase in apparent frequency toward 
present day. However, this does not occur until SD values exceed 
6.8 x 107°, a value eight times greater than that measured in the 
study area for landslides ~1700 CE (fig. $2). 

Simulated landslide chronologies are run for the past 5000 years 
and consist of two components: a background rate of sliding and 
pulses of coseismic landslides. Although we only examine landslide 
frequency over the past 1000 years, we extend the simulated chronology 
to 5000 ybp to eliminate any effects from old landslides with errone- 
ously young estimated ages. To build these chronologies, we first 
generate 10 coseismic landslide pulses at the times of past widespread 
turbidite events T1 (1700 CE) to T10. These widespread turbidite 
triggering events, identified by Goldfinger et al. (5), are thought to 
correspond to past, full-length CSZ ruptures with estimated magni- 
tudes of 8.7 to 9.1. The number of landslides triggered during each 
earthquake can be adjusted, allowing for any proportion of coseismic 
landsliding to be examined. Next, we prescribe an adjusted back- 
ground rate of landsliding, in slides/year. This represents a spatial 
and temporal average across the entire study area of the rate of non- 
seismic landslide occurrence per year. At each time step in the 
chronology, a randomly generated, normally distributed set of land- 
slide roughness values corresponding to the time step age is added. 
The background rate is constant for each model run but is scaled 
depending on the number of coseismic landslides triggered in earth- 
quake pulses so that, after adjusting for preservation bias, the total 
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number of simulated landslides estimated to have ages of less than 
1000 ybp is equal to the total number observed: ~2700 landslides. 

Because of the error associated with converting a roughness value 
to an age, coseismic landslide pulses have tails that extend past 1000 ybp, 
so the prescribed total number of coseismic landslides triggered by 
the three large earthquakes in the past 1000 years is greater than the 
total number estimated to be less than 1000 ybp. To directly com- 
pare simulated landslide chronologies to our observed chronology, 
we normalize for the same number of landslides with estimated ages 
less than 1000 years, even though this resulted in a variable number 
of prescribed landslides for each simulation. For example, for the 
simulated 50% coseismic landsliding scenario, 600 coseismic land- 
slides are prescribed per earthquake for a total of 1800 coseismic 
landslides in the past 1000 years, while only 1317 are estimated to 
have ages less than 1000 ybp (Fig. 6B). To account for the preserva- 
tion bias in the observed frequency, the total number of landslides 
for each simulated chronology is scaled as a function of time. We 
adopt a normalized exponential fit to the observed frequency data 
as this scaling function, which allows for a direct comparison be- 
tween simulated and observed frequency plots. 

We show results from simulations of 0% coseismic triggering, 
50% coseismic triggering, and 100% coseismic triggering, correspond- 
ing to 0, 600, and 1300 coseismic landslides triggered per CSZ earth- 
quake, respectively. Each of these three sets of conditions is modeled 
10* times, and then the ages are calculated for each run using the 
same best-fit roughness-age regression used to calculate observed 
ages. Last, the mean frequency as well as 5th and 95th bounding 
percentiles for all model runs are calculated and plotted (Fig. 6). 


Statistical analysis of age estimates and sources 

of uncertainty 

The total error in our estimated landslide ages stems from natural 
variability in initial roughness conditions, uncertainty in roughness 
calculations, uncertainty in interpreting radiocarbon-based landslide 
ages, and the inherently subjective nature of mapping landslides by 
eye. We attribute most of the uncertainty in our surface roughness 
calculations to variations in ground point density from the lidar point 
clouds used to generate the digital elevation models used in this study, 
an effect that has been noted in previous work (46). To quantify this 
effect and assign error bars to points on the roughness-age curve, 
we compare the roughness of two lidar datasets taken in 2011 and 
2015 of the same region within the OCR (Supplementary Materials). 
We assign the +26 m’' of the observed error between repeat lidar 
datasets (9.6 x 107+) as the +26 m™! for each point used in our 
roughness-age regression and plot this value as horizontal error bars 
on landslides of known age (Fig. 3). 

Because our age estimations are determined by the parameters 
of the best-fit regression, we needed to account for uncertainty in 
the regression itself. Uncertainty in the regression would primarily 
arise from uncertainty in the roughness-age data points used in the 
fitting process. To address the effect of error in these, we use a boot- 
strap analysis. For 10* runs, the roughness and age for each of the 
14 landslides of known age were randomly selected at a frequency 
dependent on the respective uncertainty distribution for each land- 
slide, and a best-fit exponential regression was fit. Each bootstrapped 
regression is then used to estimate landslide ages and generate land- 
slide frequency versus time plots for both the measured roughness 
values in the study area and the roughness values generated in our 
simulated coseismic landslide datasets (Fig. 5). This allows us to in- 
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terpret trends in our observed frequency in a more statistically 
robust manner. 


SUPPLEMENTARY MATERIALS 


Supplementary material for this article is available at http://advances.sciencemag.org/cgi/ 
content/full/6/38/eaba6790/DC1 
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